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ABSTRACT 


V  r — v  ^  \ 

In  a  recent  paper,  Y.  Renardy  and  D.  D.  Joseph  study  the  Benard  problem 
for  two  layers  of  different  fluids  lying  on  top  of  each  other  and  bounded  by 
walls.  Their  study  shows  that.  In  contrast  to  the  Benard  problem  for  one  fluid, 
the  onset  of  Instability  can  be  oscillatory.  The  number  of  parameters  Involved 
In  the  problem  Is  large,  and  there  Is  yet  no  comprehensive  picture  of  when  the 
Instability  Is  oscillatory  and  when  It  Is  not.  The  study  of  limiting  cases, 
accessible  by  perturbation  methods,  may  be  helpful  In  this  respect.  In  this 
-paper,  an  analysts  Is  given  for  the  case  when  the  properties  of  the  two  fluids 


are  nearly  equal  and  the  fit  ids  are  allotted  to  slip  at  the  boundaries. 
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SIGNIFICANCE  AND  EXPLANATION 


Flows  involving  two  incompressible  viscous  fluids  exhibit  nonuniqueness  in  the  sense 
that  many  interface  positions  are  allowed  when  their  densities  are  equal.  Two-fluid  flows 
also  have  quite  different  dynamical  features  from  one-fluid  flows.  The  one-fluid  Benard 
problem  in  which  the  fluid,  lying  between  parallel  horizontal  plates,  is  heated  from  below 
has  a  static  solution  for  which  a  linear  stability  analysis  yields  no  complex  eigenvalues.  On 
the  other  hand,  the  two-fluid  problem  yields  complex  eigenvalues.  In  this  paper,  we  use 
perturbation  methods  to  examine  the  conditions  under  which  such  time-periodic  instabil¬ 
ities  occur.  This  may  have  application  to  the  theory  of  convection  in  the  Earth’s  mantle, 
which  is  sometimes  based  on  the  assumption  that  convection  takes  place  in  chemically 
uniform  layers. 


The  responsibility  for  the  wording  and  views  expressed  In  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  of  this  report. 


PERTURBATION  OF  A  MULTIPLE  EIGENVALUE 
IN  THE  B&NARD  PROBLEM  FOR  TWO  FLUID  LAYERS 

Yuriko  Renardy  and  Michael  Renardy 


1.  INTRODUCTION 

In  the  Blnard  problem  for  one  fluid,  “exchange  of  stabilities’'  holds  for  a  variety  of 
boundary  conditions1,  whether  the  fluid  is  bounded  by  walls  or  by  stress-free  surfaces, 
or  by  a  wall  below  and  by  a  gas  above.  Critical  eigenvalues  are  not  oscillatory  in  time. 
The  consideration  of  additional  effects  can,  however,  introduce  “overstability”,  that  is, 
critical  eigenvalues  that  are  oscillatory  in  time.  Overstability  occurs,  for  example,  if  there 
are  temperature-dependent  surface  tension  gradients3,  if  there  is  a  temperature-dependent 
solute  gradient3,  in  mixtures  of  superfluids4,  or  if  two  fluids  lie  in  layers  with  different 
solute  gradients**6.  On  ther  other  hand,  the  Bdnard  problem  with  two  layers  of  fluids  with 
different  thermal  and  mechanical  properties,  without  surface  tension  gradients  or  solutes, 
has  only  recently  been  examined  for  the  possibility  of  overstability7.  Moreover,  Busse  has 
suggested  that  such  a  model  may  explain  certain  features  of  mantle  convection,  such  as 
the  sise  and  aspect  ratio  of  the  convection  cells  8  and  the  possibility  of  time-periodic  9 
flows. 

In  a  recent  study,  Renardy  and  Joseph  7  show  that,  in  contrast  to  Etenard  convection 
for  one  fluid,  the  linear  stability  problem  for  two  fluids  in  layers  is  not  self-adjoint,  so  that 
complex  eigenvalues  are  possible.  They  then  solve  the  eigenvalue  problem  numerically  and 
find  an  oscillatory  onset  of  instability  in  a  situation  where  the  two  fluids  are  only  slightly 

Sponsored  by  the  United  States  Army  under  Contract  No.  DA  AG  29-80-C-0041  and 
supported  in  part  by  the  Centre  for  Mathematical  Analysis  at  The  Australian  National 
University. 


different.  It  would  be  interesting  to  know  what  determines  whether  the  onset  of  Blnard 
convection  is  stationary  or  oscillatory.  However,  an  attempt  to  answer  this  question  by 
numerical  calculations  is  not  feasible  because  of  the  large  number  of  parameters  involved. 
We  think  that  it  is  therefore  helpful  to  first  find  out  what  happens  in  certain  limiting 
situations  accessible  by  perturbation  methods. 

In  this  paper,  we  look  at  two  horizontal  fluid  layers,  lying  between  parallel  boundaries 
which  are  kept  at  constant  but  different  temperatures.  In  order  that  the  solution  to  the 
unperturbed  problem  be  available  in  closed  form, we  require  that  slip  boundary  conditions 
apply,  that  is,  the  normal  velocity  and  shear  stress  vanish.  When  the  properties  of  the 
two  fluids  are  the  same,  the  eigenfunctions  at  criticality  and  the  critical  Rayleigh  number 
can  be  determined  exactly  1  and  there  is  a  double  eigenvalue.  When  the  properties  of 
the  fluids  are  slightly  different,  we  obtain  a  regular  perturbation  expansion  for  the  double 
eigenvalue,  for  which  we  calculate  the  two  leading  terms.  This  allows  us  to  investigate  the 
following  two  questions: 

1.  Which  perturbations  lead  to  eigenvalues  with  nonvanishing  imaginary  parts? 


2.  Which  perturbations  stabilize  and  which  destabilize  the  flow? 


2.  GOVERNING  EQUATIONS 


We  consider  a  fluid,  filling  the  space  between  two  parallel  boundaries  of  infinite  extent 
in  the  (x*,z*)-plane.  Asterisks  denote  dimensional  variables.  The  upper  boundary  at 
z'  =  /'  is  kept  at  a  constant  temperature  T^,  and  the  lower  plate  is  kept  at  a  higher 
constant  temperature  Tq  +  A  T*.  At  temperature  T^,  the  fluid  has  a  coefficient  of  cubical 
expansion  &  ,  thermal  diffusivity  k,  thermal  conductivity  k,  viscosity  /*,  kinematic  viscosity 
1/  and  density  p.  We  use  the  Boussinesq  approximation  in  the  Navier-Stokes  equations, 
that  is,  the  density  in  the  buoyancy  term  is  approximated  by  p(l  -  &{T*  —  Tq)). 

Following  Drazin  and  Reid  1,  we  introduce  dimensionless  variables  (without  asterisks) 
as  follows: 

(x,z)  =  (x',z*)/r, 

t  =  Kf/V\  (1) 

u  =  u ’/'/«, 

T  =  T*/AT\ 
r=,TinpK*). 

Here,  u’  =  (u",w')  is  the  velocity,  p"  the  pressure  and  T"  the  temperature.  We  define 
a  Rayleigh  number 

R  =  ga£ari'zl{Kv), 

where  g  denotes  gravitational  acceleration,  and  a  Prandtl  number 


We  study  the  linear  stability  of  the  static  solution 

u  —  0,  T  =  T0  +  1  -  z. 

If  disturbances  are  proportional  to  exp(at  +  tax),  then  the  following  eigenvalue  problem 
arises  for  the  velocity  u  and  perturbation  6  to  the  temperature: 

<70  =  A0  +  w, 

au  =  -  V  P  +  RP&e.z  +  P  A  u  ,  (2) 

V  •  u  =  0. 

At  the  boundaries,  w  =  0  and  0=0.  At  solid  boundaries,  we  would,  in  addition,  have 
u  =  0.  However,  we  prefer  to  look  at  the  case  of  slip  boundary  conditions,  where  it  is 
required  that  the  shear  stress  be  zero,  or  equivalently,  du/dz  =  0.  Although  this  is  not 
physically  realistic,  this  has  the  advantage  that  the  eigenvalue  and  eigenfunction  at  the 
critical  Rayleigh  number  are  known  in  closed  form.  It  is  well  known  1  that  all  eigenvalues 
are  real,  and  the  critical  case  o=0  occurs  first  at  R  =  27tt4/4.  The  wavenumber  of  the 
critical  mode  in  the  x-direction  is  a  =  2-1/2jt. 

We  perturb  this  marginally  stable  eigenvalue  by  the  introduction  of  a  second  fluid, 
whose  properties  are  slightly  different.  At  rest,  this  introduces  an  interface  at  z  =  /)  as 
indicated  in  Fig.  1.  Subscripts  1  or  2  on  the  physical  parameters  will  now  denote  fluid  1 
or  2,  respectively.  There  are  6  dimensionless  ratios: 


Sketch  of  B4nard  problem  perturbed  by  a  second  fluid 
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The  non-dimensionalization  of  the  equations  and  the  definitions  of  Rayleigh  and  Prandtl 
numbers  will  be  based  on  fluid  1;  for  example,  R  =  g<k\AT*l*3 I[kiV\).  We  include  a  sur¬ 
face  tension  between  the  fluids,  described  by  a  dimensionless  parameter  S  —  Sml‘ 
where  5 *  is  the  dimensional  surface  tension  coefficient. 

The  linearized  eigenvalue  problem  reads  now  as  follows: 

For  0  <  z  <  li. 


<70  =  wAi  +  AO, 


au  =  -dp/dx  +  P  A  u, 


aw  =  —dp/dz  +  RPQ  +  P  A  w , 


du  dw 

di  +  fa  =  ’ 


where  At  =  For  /i  <  z  <  1, 


(4) 


<70  =  Wy42  +  —  A  0, 

7 

dp  r 

<7tt  =  -T-~  +  -P  AU, 

ox  m 

dp  RP^  r  „ 

aw  =  -r  — ^  +  — t-0  +  —  P  Aw, 
az  p  m 


(5) 


v  •  •  v  ■,*  V  V  V  /  * 


where  Aj  =  (At. 


du  dtu 

di  +  HI  ~  ’ 


The  boundary  conditions  at  z= 0  and  1  are  as  before: 


du 

0  =  w  =  — -  =  0. 

dz 


(6) 


The  interface  is  perturbed  to  the  position  z  =  ly  +  h(x,t).  The  interface  conditions 
linearized  about  z—ly  are: 

continuity  of  temperature:  [  [9  ]  ]  =  h  [  [A  ]  ] , 
continuity  of  heat  flux:  [  [kff  ]]  =0, 

continuity  of  velocity:  [[«/]]  =  [  [u  ]  ]  =0,  (7) 

continuity  of  shear  stress:  [[p(f*j  +  )  ]]=0, 

balance  of  normal  stress: 


dwi  1  du>2  „  d^h 


Here,  [[.  ]]  denotes  the  jump  of  a  quantity  across  the  interface,  for  example, 


and  we  have  set 


The  kinematic  free  surface  condition  yields,  at  z  =  1 y. 


We  now  proceed  as  follows.  We  fix  R  =“^*4  and  keep  P  and  l\  fixed  but  arbitrary. 
We  consider  only  solutions  to  (4)  -  (8)  proportional  to  exp(taz)  with  a  fixed  at  2 -1/,2x. 
We  introduce  a  small  parameter  e  and  regard  1  -  m,  1  -  r,  1  -  7, 1  —  f,  1  -  0,M\  and  S  as 
small  quantities  proportional  to  c;  that  is,  we  set 

1  -  m  =  me, 

1  -  r  =  f£, 

1  — y  = 

1  -  ?  = 
l-0  =  fc, 

YIk*  P  f 

Mi  =  Mjf,  A/j  =  — - — (-  +  Mf  +  ^))» 

5  =  Se. 

At  e  =  0,  there  is  an  algebraically  two-fold,  but  geometrically  simple  eigenvalue  a  =  0. 
One  eigenvalue  arises  from  the  first  criticality  of  the  Benard  problem  for  one  fluid  and  the 
other  arises  from  the  presence  of  the  interface.  For  small  £,  this  eigenvalue  is  perturbed 
into  two  eigenvalues  which  can  be  expanded  in  powers  of  e l^2.  The  purpose  of  the  following 
analysis  is  to  find  the  coefficients  of  e and  £  in  this  expansion. 


8 


3.  PERTURBATION  OF  MULTIPLE  EIGENVALUES 


The  perturbation  expansion  for  a  double  eigenvalue  is  not  a  simple  series  expansion 
in  (.  In  addition,  eigenfunctions  need  not  have  series  expansions  in  powers  of  t.  The 
perturbation  of  multiple  eigenvalues  is  discussed  in  Ref.  10,  chapter  IV  ,  §1.  The  procedure 
involves  the  generalized  eigenspaces  belonging  to  the  eigenvalue  o  =  0  for  the  unperturbed 
problem  and  its  adjoint,  and  does  not  require  finding  the  eigenspaces  of  the  perturbed  O(e) 
problem  at  all. 

We  now  quote  the  pertinent  results  from  Ref.  10.  Suppose  o0  is  an  algebraically 
2-fold  eigenvalue  of  a  matrix  Lq.  Let  {01,02}  be  a  basis  for  the  generalized  eigenspace 
of  Lq  with  eigenvalue  00,  and  let  {61,62}  be  a  basis  for  the  generalized  eigenspace  of  Lq 
(the  adjoint  of  L0)  with  eigenvalue  do  (the  overbar  here  denotes  the  complex  conjugate). 
Let  Lq  be  perturbed  into  L(c)  =  Lo  +  <Li(e)  +  0(c2)  with  L\  depending  smoothly  on  c. 
Then  the  perturbed  eigenvalues  o  are  given  by  the  zeros  of  the  determinant  of  a  matrix 
¥t;(e, o), i,j,  =  1,2,  which  represents  to  O(e)  the  projection  of  L(e)  -  o,  first  onto  the 
eigenspace  of  the  unperturbed  problem  and  then  onto  the  adjoint  eigenspace: 

*o(£»a)  =  (6,,  {Lo  +  tLi(e)  -  o)aj )  +  0(c2).  (9) 

Some  care  must  be  taken  when  this  result  is  applied  to  unbounded  operators  in  infinite 
dimensional  spaces,  for  example,  differential  operators.  Such  an  operator  has  a  “ domain ” 
that  is  specified  not  only  by  smoothness  requirements  on  the  function  but  also  by  the 
boundary  conditions.  If  the  domain  of  the  operator  that  is  being  perturbed  depends  on  e, 
we  cannot  apply  (9);  the  domains  of  L(<)  and  Lo  may  be  different,  and  their  combination 
would  not  make  sense.  We  can,  however,  circumvent  this  problem  by  not  looking  at  the 


differential  operator  itself,  but  at  its  resolvent  (L(e)  -  A/)' 1  where  A  is  not  an  eigenvalue  of 
L(f).  The  domain  of  this  does  not  depend  on  c  and  we  will  need  to  redefine  accordingly. 
As  will  be  seen  from  the  following,  the  resolvent  itself  does  not  ever  need  to  be  computed. 

In  order  to  make  these  ideas  more  precise,  we  introduce  some  notation.  Let  X  denote 
the  set  of  functions  (0,u,tu,/i).  We  introduce  an  inner  product  by 

i.2»/a 

(Xl,X7)  =  /  /  ©102  +  U1U2  +  W1W2  dzdx 

JO  Jz= 0 

[2* /a  1 1 

+  /  /  ©i©2  4-  ®itt2  +  W\U)2  dzdx 

Jo  Jg=h 

[2r/a _ 

+  /  h\h,2  dx  (10) 

Jo 

to  generate  a  Hilbert  space.  In  this  Hilbert  space,  we  consider  the  subspace  determined  by 
the  “Hodge  projection”  (see  space  H  in  Theorem  1.4,  Ref.  11),  that  is,  by  the  conditions 
that  the  velocity  field  be  divergence-free,  that  the  vertical  velocity  vanish  at  the  bound¬ 
aries,  and  be  continuous  across  the  interface.  By  L(e)X  we  denote  the  right  hand  sides 
of  equations  (4), (5)  and  (8).  We  regard  L(e)  as  an  operator  in  the  subspace  so  that  the 
conditions  on  w  in. (6)  and  (7)  and  the  normal  stress  balance  in  (7)  are  an  integral  part 
of  the  definition  of  L(e).  The  domain  of  definition  of  L(f)  is  determined  by  the  rest  of  the 
boundary  conditions  in  (6)  and  (7),  which  we  write  in  the  form  B(e)X=0.  The  range  of 
the  operator  L(c)  must  satisfy  the  following  conditions  in  order  for  the  pressure  p  occuring 
on  the  right  sides  of  (4)  and  (5)  to  be  determined  as  a  function  of  X:  The  “velocity  part” 
of  L(()X  must  be  divergence  free,  the  vertical  velocity  must  vanish  on  the  walls  and  be 
continuous  across  the  interface,  and  the  jump  in  p  across  the  interface  must  be  given  by  the 
normal  stress  balance.  Thus,  the  problem  we  wish  to  solve  is:  for  small  c,  find  a  satisfying 


»  •* »  * t  **»•*.  *" «  ^ .  **.  *'»  •*•  •*.  *"• 
*  .*•  .'  .'  /*  !  •  -  >  >  .*•  /. 
syvV/VvV.-vv‘  >  • 


l  /ss1 


'.s  .■•v-y-vv ,  '•].**] 
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V  W'  , 


V 

V' 

Ov 


L(e)X  =  aXy 


B(e)X  =  0, 


L(t)  =  Lq  +  iL\  +  0(c2), 


Y'.- 


I 


Explicitly, 


in  fluids  1  and  2,  and 


in  fluid  1  and 


fl(f)  —  Bo  +  (Bj  +  0(e2). 


LnX  = 


A©  +  w  \ 
-||  +  PA» 

-§ £  +  RPe  +  PAw 
w 


LtX  = 


_2e 

gt 

-§£ 
da 

\  o  J 


^A0-/ifu> 

'If  +  (*.-*')/’  A«-H 

f  §*  +  $RPB  +  (m  -  f)P  A  ti;  -  || 
0 


in  fluid  2,  where  p  denotes  the  (7(c) -perturbation  to  the  pressure, 


B0X  = 


(  0i  -  02  at  z  =  l\  \ 
*t  * = >• 

ui  —  U2  at  z  =  l\ 

§JLl  _  2*2  at  2  =  1. 
da  da  M  2  *> 

0  at  z  =  0, 1 
§*  at  z  =  0, 1 


J 


and 


B\X  = 


-hf  at  z  =  /i 
-f  ^  at  z  =  /, 

0  at  z  =  /j 
+  fe)at2  =  (, 
0  at  z  =  0, 1 
0  at  z  =  0, 1 
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With  the  above  definitions,  we  are  now  ready  to  look  at  the  resolvent  of  L(c)  and  then 
to  redefine  the  matrix  in  (9).  Since  c  is  small,  the  eigenvalues  of  L(e)  are  close  to  those  of 
Lq,  so  that  the  A  in  the  resolvent  should  be  chosen  well  away  from  sero.  We  choose  A  =  1. 
Hence,  instead  of  looking  at  (11)  directly,  wc  study  the  equivalent  problem 


(1(c)  -  1)-*X  =  (a  -  l)'1*  =:  dX 


and  perturb  around  a  —  -1.  We  note  that  the  definition  of  (L(c)  -  l)-1  already  in¬ 
corporates  the  boundary  conditions.  The  relation  (9)  is  applied  to  this  problem.  The 
determinant  of  the  matrix 


*«(«»*)  =  (fit,  ((£(«)  -  1)  1  -  o)aj)  +  0(c2),  ij  =  1,2, 


•  ■  ;1 
>  \-  ■  :■ 


where  the  6,  and  ay  are  as  before,  is  set  equal  to  zero.  We  will  require  an  expansion  of 
the  resolvent  in  powers  of  c  in  order  to  carry  out  the  calculation  of  4\y.  We  note  that  the 
inverse  of  L'o  ~  1  is  defined.  Since  parts  of  the  calculation  are  rather  lengthy,  they  will  be 
organized  into  several  appendices  of  this  paper. 

We  first  have  to  find  the  boundary  value  problem  adjoint  to  (11).  This  is  done  in 
Appendix  A.  Then  we  determine  the  generalized  eigenvectors  at  e=0  for  both  (11)  and  the 
adjoint  problem.  These  are  denoted  by  a,,  02  and  61,62,  respectively,  and  satisfy: 


*  .  - 


loo,  —  0,  B0a,  =  0, 


LqC2  =  Oj,  B 0O2  =  0, 


Lobi  =  0,  B^b  1  =  0, 


Z-0^2  =  Bqbi  =  0. 


* ,  *  •  *  1  •  •,*  •,  /,  *.  *«  /•  ,■*,  * .  •  .* ,» ,*•  ,*• .  • .  -  .*•  /■ 


*  •  *  • « •  %  *  .««•.* 
r,»  •  •  ■ 


V-^ 
v  v*. 
VV 
..... 


V- 


The  generalised  eigenvectors  are  determined  in  Appendix  B.  In  order  to  apply  formula 


(12),  we  must  determine  the  expressions 

(MLW-l)-1.,-)  04) 

to  first  order  in  t.  To  facilitate  this  calculation,  we  introduce  x°  and  xj  defined  by 

{L(t)  -  1  )_,Oy  =  xj  +  IX j  +  0(|J).  (15) 

Equating  the  coefficients  of  equal  powers  of  c,  we  find  the  equations  governing  x°  and  xj 

(Lq  -  l)xj  =  Oy, 

Box?  =  0,  (16) 

and 

Z»lXy  +  (Lo  -  l)xy  =  0, 

fl.xj  +  B0xj  =  0.  (17) 

From  (16),  we  immediately  find  x°  =  -«!,  x§  =  -oj-oj.  We  will  not  need  the  solutions  xj 
to  the  perturbation  problem  (17)  but  only  certain  inner  products  involving  them,  namely 
(6,,Xy).  This  is  seen  from  (12)  and  (15): 

♦<>(*»*)  =  <fc,*y)  +  <<*«,*j)  -  o(bi,aj)  +  0{(2).  (18) 

We  calculate  (6,-,xj)  from  (17): 

(5|,Xy)  =  (6,,L,Xy)  +  (6i,L0Xy},  (19) 
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v.v 


V 
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and  an  integration  by  parts: 


{bi,Loxj)  =  (Lo6„xj)  +  boundary  integrals ,  (20) 

where  the  boundary  integrals  are  evaluated  using  the  second  part  of  (17).  (The  boundary 
integrals  would  vanish  if  B()Xj  were  zero.)  Details  of  these  calculations  are  in  Appendix  C. 

In  the  following  section,  we  discuss  the  solution  of  the  eigenvalue  perturbation  prob¬ 
lem,  which  now  reads  detftij  =  0,  where  ♦  is  given  by  (18). 
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4.  RESULTS  AND  DISCUSSION 


We  set  d  —  -1  +  T\f}! 2  +  r2e  +  0(e3/2).  Since  a  —  1  +  1/d, 

a  =  -rlv/t  -  (t,2  +  t2)(  +  0(f3/2). 

For  the  first-order  term,  we  find  (see  equations  C2,  C5,  C7  and  CIO  of  Appendix  C) 

r?  =  -  %r+  ’)'  (»*?' +  J **l-  +  W  +  f)|  +  s>-  (2,) 

The  sign  of  this  quantity  determines  whether  the  eigenvalue  splits  into  real  values  or 
complex  conjugates.  We  note  that,  in  the  situation  studied  by  Renardy  and  Joseph  7,  f, 
f  and  S  were  0,  and  $  was  positive.  (Their  boundary  conditions  are  different  from  ours, 
nevertheless,  we  expect  some  qualitative  similarity).  They  find  complex  eigenvalues,  in 
agreement  with  the  prediction  of  (21).  It  is  also  worth  noting  that  (21)  does  not  involve 
the  viscosity  ratio,  and  that  the  sign  does  not  depend  on  the  Prandtl  number. 

If  r?  is  positive,  we  always  have  instability.  However,  if  r2  is  negative,  we  have  to  find 
r2  +  r2  in  order  to  determine  whether  the  perturbations  are  stabilizing  or  destabilizing. 
Unfortunately,  the  formula  for  rj  is  not  nearly  as  simple  as  that  for  T\  (see  equations  C5, 
C6,  Cll,  and  C15  of  Appendix  C).  We  find 

r?  +  r»  =  +  -*ft)  +  (A3Ti+  yP5)sinx/, 

+  t(1  “  HI*)fcos#l|  +  Pf-^ReaI[c,(Q?  -  ft2l2)iQx coshQi/i) 

s  ft  * 

+(1  +  n2(C^~T  +fi*in#/i) 

s  sin  fti  i 

+(JWT,  +  ~P§)(2Real\ci  sinhQ|/i|  +  (1  +  4)^  cos*/,) 


(22) 


+^(‘  ■  3 MiT?)(*  +  (L7FL('":o‘'',+  £»)• 

where  Mj  is  defined  after  equation  (8),  ct  and  Q i  are  defined  in  Appendix  B. 


We  examine  some  limiting  cases  below.  If  the  Boussinesq  approximation  is  to  be 
justified,  then  diA T*  must  be  small1,  e.g.  a  =  5.10 ~4K~l  and  AT*  <  10K.  For  the 
purpose  of  computing  graphs,  we  have  taken  diAT*  =  0.001.  As  P  — ►  oo,  both  r,2  and 
t\  +  r2  approach  finite  limits. 

Case  i.  <S^0.f=f  =  fl  =  rh  =  0.  If  S  <  0,  then  r2  >  0  so  instabilities  are  exponen¬ 
tial  in  time  and  leads  to  mixing.  If  5  >  0,  then  r2  <  0  so  the  surface  tension  opposes  the 
tendency  for  convection.  However,  whether  this  leads  to  a  growth  or  not  depends  on  the 
sign  of  (22)  which  depends  nonlinearly  on  /i  and  the  Prandtl  number.  Fig.2  shows  this 
dependence  to  be  r2  +  r2  >  0  so  that  a  is  stable  and  oscillatory. 

Case  ii.  f  ^  0.  S  =  (  =  3  =  th  =  0.  Since  djA T'  is  small,  r2  and  f  have  the  same 
signs:  if  Fluid  1  is  the  less  dense  fluid,  a  convective  instability  takes  place;  if  Fluid  2  is  the 
more  dense  fluid,  stable  wavy  modes  occur.  The  graph  of  r2  +  r2  versus  l\  for  f  <  0  is 
similar  to  Fig.2. 

Case  iii.  0^Q.S  —  =  0.  The  signs  of  r2  and  $  are  opposite.  A  heuristic 

explanation  is  as  follows.  In  the  Boussinesq  approximation,  the  density  at  any  depth  / 1  in 
the  buoyancy  term  of  the  momentum  equations  is  approximated  by  p,(l  -  d,  (T*  —  Tq  )),  i  = 
1,2  where  (Tm  -  Tq)/  AT'  =  l-?  +  0.  Hence,  at  a  depth  z  which  is  close  to  the  interface 
z  —  / 1 ,  the  densities  are  approximately  Pi(l  -  d,/2).  Therefore,  if  <  0,  the  lower  fluid 
is  less  dense  than  the  upper  fluid,  locally  at  the  interface,  so  that  r?  >  0.  On  the  other 


0.4  0.6  , 

Depth  of  lower  fluid  c, 


hand,  if  $  >  0,  then  the  lower  fluid  ie  the  more  dense  and  rf  <  0.  Fig.3  presents  r?  +  r2 
versus  fj  for  4  =  showing  that  for  a  wide  range  of  ii  close  to  0  and  of  Prandtl  numbers 
greater  than  0,  there  are  unstable  oscillatory  modes.  These  modes  are  Hopf  bifurcations. 


Case  iv.  (^0,S  =  f  =  3  =  fh  =  0.  The  dependence  of  to/j  is  through  -f  sin  2x7  §. 
Hence,  if  the  thicker  layer  has  the  lesser  coefficient  of  conductivity,  (i.e.  if  f  >  0  and 
0.5  <  / 1  <  1;  or  if  f  <  0  and  0  <  l\  <  0.5)  then  convective  instability  results.  If  the 
thicker  layer  has  the  greater  coefficient  of  conductivity,  then  time-periodic  modes  occur, 
but  whether  these  would  be  stable  or  not  depends  on  the  sign  of  rf  +  r2.  Fig.4  shows 
a  wide  range  of  unstable  oscillatory  modes  for  0  <  l\  <  0.5  when  f  =  1.  This  graph  is 
antisymmetric  about  7j  =  0.5,  as  can  be  deduced  non-trivially  from  the  equations.  We 
note  that  when  the  fluid  with  the  greater  conductivity  occupies  most  of  the  flowfield,  the 
stability  of  the  time-periodic  eigenvalue  a  depends  on  the  Prandtl  number:  if  the  Prandtl 
number  is  very  small,  the  oscillatory  modes  are  stable,  and  if  the  Prandtl  number  is  well 
away  from  0,  the  oscillatory  modes  are  unstable. 

Case  v.  tn  ^  Q.(  —  S  =  f  =  5  =  0.  We  find  that  rf  =  0,  and  Fig.5  shows  r2  +  r2 
versus  /j  for  various  Prandtl  numbers.  If  Fluid  1  is  the  less  viscous  fluid,  then  there 
is  stability,  and  if  Fluid  1  is  the  more  viscous  fluid  then  there  is  convective  instability, 
regardless  of  the  depth  ratio,  in  agreement  with  heuristic  expectation. 


=  f  =  S  =  f  =  3  =  0.  The  ratio  7  of  thermal  diffusivities  plays  a 


similar  role  here  to  the  ratio  m  of  viscosities,  which  measures  the  rates  of  diffusion  of 
momenta.  As  in  Case  v  above,  rf  =  0  and  the  graph  of  r 2  -f  r2  versus  /]  is  shown  in  Fig.6. 
If  Fluid  1  has  the  lesser  coefficient  of  thermal  diffusivity,  then  a  is  stable.  If  Fluid  1  has 
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Graph  of  r?  +  r2  versus  l\ for  various  Prandt!  numbers;  /f  =  l,  5  =  f  =  f  =  if»*^  =  0, 


iki  • 


Tl"2+T2 


Tl"2+T2 


the  greater  coefficient  of  thermal  diffusivity,  then  convective  instability  ensues,  regardless 
of  the  depth  ratio,  as  expected. 
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APPENDIX  A:  The  adjoint  problem  for  f =0 


We  denote  the  domains  occupied  by  the  two  fluids  by 


Hi  =  {0  <  x  <  2*/a,0  <z<l\ } 


and 


fl2  =  {0  <  x  <  2ir/a,l\  <  z  <  l}. 


We  denote  the  interface  by  I,  and  the  lower  and  upper  boundaries  by  T j  and  r2,  respec¬ 
tively.  Let  Xi  =  (0,  u,w,h)  and  X2  =  (G*,u’  ,w' ,  h‘).  Asterisks  denote  the  adjoint.  We 
have 


{X2,L0Xi)=  f  0*  A  8  +  0*tu  -u*^~  +  Pu'  A  u  -  te*  ^  +  RPw’Q  +  Pw*  A 
J  n,  ox  dz 


w 


+ 


/  0*  A  0  +  0*u>  -  u* ^  +  Pu*  A  u  -  ti)‘ ^  +  PPtB'0  +  Pib*  Aid 
Jn-, 


dx 


dz 


+  J  h*W\.  (/li) 

We  integrate  by  parts  and  obtain,  using  the  divergence  condition  div  u  =  div  u"  =  0, 


(X2,L0X i>=  f  0(A8'  +  RPw*)  +  u(P  Afi"  -  ~z~)  +  W{P  A  ii&"  —  +  ©') 

J  n,  ox  dz 


+ 


[  0(A0‘  +  RPw*)  +  u(P  Au‘-  +  u/(P  A®'  -  +  0-) 

dx 


dx 


+ 


/ 

xr, 


d©  d©‘  ,du  du*  tdw  dtD*  .  „ 

-_  +  e___pu  ^  +  Pu— -Ptf  _  +  Pu;__  +  tDp_u,p 


+ 


+  JiK.a,  +  [[e-£-e^  +  P.-*“ 


**•1 

dz  dz 


dw  dw *  , , 

+P6)'  -  *>+ 


M2) 


From  this,  we  immediately  read  off  the  adjoint  differential  operator  to  be  given  by 


l;,x2  = 


A  e+RPw' 
p  A«--& 
p Aw-^f  +  e*  ' 


Moreover,  since  Xj  satisfies  the  boundary  conditions  (6)  on  f  i  and  r2,  the  integrals  over 


these  boundaries  vanish  if 


du’ 

9‘  =  w'  =  ~  =o. 
dz 


Into  the  interface  term  in  (A2),  we  add 


_  A.du  dw .  .di*  dw\ 

-p»(*  +  JI)  +  p"(te+ar> 


which  is  zero.  We  integrate  the  x-derivative  by  parts  and  use  periodicity.  This  yields 
ft.  „  ./d«  n  ,dfi *  diD\ 

+  [[•«;-  e«T  +  <51 +  -  Pu(  a?  +  *r> 


From  this  we  find  the  adjoint  interface  conditions: 


[[»']]  =0, 

i&O’ 

[[•■]]  =0. 

[Ml  =0. 

rr3u'  dtlT,, 

!1-5T+  *rll  =  0> 

p  -2pd-£rW  +*'  =  «• 


> * •. ^ *, 


APPENDIX  B:  Eigenfunctions  of  the  unperturbed  problem 

If  e=0,  the  variable  h  does  not  occur  in  the  right  hand  sides  of  (4),  (5)  ,(8),  or  in  the 
interface  conditions,  and  we  have  the  eigenfunction 


„  _  -IOI 

a\  -  t 


The  equations  (4)-(7)  are  precisely  those  characterizing  the  one-fluid  Benard  problem.  The 
eigenfunction  for  this  problem  1  now  yields  a  generalized  eigenfunction: 


/  *» 

r3+a3 

Mxaz  tv 

t  “  cos  irz 

a  2  =  -  “ 

sin  irli  sin?rz 


The  adjoint  equations  agree  with  the  one-fluid  Bdnard  problem  if  we  set  h'  =  0.  Thus  the 
eigenfunction  of  the  one-fluid  Bdnard  problem  yields  the  eigenfunction  for  the  adjoint: 


6,  =  e’ai 


|jPjr2sin  icz 

COS  7TZ 

sin  itz 


The  generalized  eigenvector  62  of  the  adjoint  satisfies 


L'0b  2  =  61, 


B'0b2  =  0. 


This  leads  to  the  equations: 


A©'  +  RPw'  =  - Pir2 e,ax  sin  it z, 
2 


=  —  e,aIcos  jtz, 
ox  a 


P  Ato'  +  0‘  -  =  e  a,sin  jtz, 

az 

du'  dw’ 

+  ~dT  ~  °* 


We  set  tt>'  =  u>(',e*a*  etc.,  and  obtain  by  combining  the  equations: 


/  ^  2\3  *  27  Q  *  9  4/-  1  \  • 

(rfz*  "  °  )  *0  +  “>o  =  ^  +  p)s«njr«. 


The  general  solution  of  this  equation  is 


Wq  =  e  i  sinhQiz  +  ejcoshQiz  +  cssinhQ2Z  +  C4CoshQjz 


1  .  1 . 

+cssin jtz  +  cocos nz  -  —  (1  +  ~)zcos7rz 

6*  r 


in  fluid  1,  and 


Wo  =  dtsinhQj(z  -  1)  +  di cosh Q  1(2  -  1)  +  d3sinhQ2(*  -  1)  +  d4 coshQ2(z  -  1) 
+d5sin *{z  -  1)  +  d&c<XK(z  -  1)  -  ^(1  +  ^Jzcosjtz  (J 


in  fluid  2,  where 


9.  =  «-W*. 


and  <t>  is  determined  by  cos  4>  =  5/\/52,  sin <j>  =  Zy/Z/52.  The  coefficients  e  1  *  cq  and 

dj  -  do  must  be  determined  such  that  the  boundary  conditions  are  satisfied.  By  using 
(B4),  we  can  show  that  the  conditions  (A4)  at  the  walls  reduce  to 


At  z=0,  this  yields 


d2  d4 

<  ~  ^  =  °- 


Cj  +  C4  +  cj  =  0, 


27 


Q l«2  +  Q\cA  -  *  C6  =  o, 


QjCj  +  Q^c  4  +  ir4ce  =  0- 


From  this,  we  obtain  -  c4  -  cq  -  0.  At  z=l,  we  find 


<*2  +  <*4  +  do  +  —(1  +  p)  =  0, 


(BIO) 


Qidj  +  Q2d4  -  jr2do  -  —  (1  +  p)  =  0, 

Q*d2  +  Q\d4  +  jr4do  +  ■g-  (1  +  p)  =  0. 


(Bll) 


This  yields  d2  =  d4  =  0,  d6  =  -  £(1  +  £).  The  first  five  of  the  conditions  (A6)  lead,  after 
eliminating  u*  and  0*  from  (B4),  to  the  conditions 

IK)]  =  [[^]]  =  [[^11  =  [[£?]] 


(B12) 


We  can  set  the  coefficient  <f5  to  zero  for  the  following  reason.  The  coefficient  rfs  multiplies 
ui  =  sin  ir ( z  -  1)  in  the  generalized  eigenvector.  The  eigenvector  6j  has  ti>  =  sin  xz.  Since 
any  multiple  of  6j,  added  to  62,  i.e.  +  C6i,  is  also  a  generalized  eigenvector,  we  choose 
C  to  be  d$.  This  essentially  gets  rid  of  d$  in  Qj,  replaces  C5  in  flj  by  c$  +  d$,  and  we 
rename  c$  4-  ds  as  C5.  We  thus  obtain  the  following  system  of  equations. 


ci  sinhQifi  +  c3sinhQ2fi  +  cssinz-fi  =  -di  sinhQi/2  -  <f3sinhQ2/2  +  decosz/2, 


CiQjcoshQi/i  +  c3(?2cosh(?2/i  +  Jrcscosz/i  =  (?i<fiCoshQ,/2 


+02*3  cosh  Q21 2  +  Trdesin  irl  2, 


m 


ciQ^sinhQj/j  +  c3Qf  sinhQ2/j  -  n ^5  sin  xlx  =  -QfdiainhQih 
-<?2rf3sinhQ2/2  -  *2dcCosxl2, 


C]Q{  sinhQi/i  +  c3Q4sinhQ2/i  +  *4c$  sinjr/j  =  -Q\dx  sinhQ|/2 

-02ds  sinh  Q2/2  +  x4d6  cos  *72,  (£13) 

cx(Q*  -  jt3Q|) cosh Q \l\  +  c3((?j  -  jr2Q®)  cosh Q  2^  +  2c5jr5coeir/| 

=  dx(Q\  -  7r2Qf)co6hQi/j  +  d3(Q2  -  *2Q2)  coshQ2/2  +  2d0*5  sin*/2. 

We  eliminate  ex  from  the  third  and  fourth  equations  by  using  the  first  equation. 

(Q*  ~  Q 1)  «nh Q*l i®3  -  (»2  -  Qi)  sin  nl ic5  +  (Q2  -  Q 2 )  sinh Q2/2  =  - (s2  +  Q\ )  cos  ir l2d6, 

(Qi(Ql  “  *2)  ~  QiQiiQx  -  »2))coshQ2/iC3  +  (2w6  -  *Q?(Qf  -  Jr2))eosjr/iC5 
+(-<?!(<?2  -  *2)  +  QaQ2(Q2  -  *2))  cosh  Q2l2d3  =  (2jt5  -  Q\*{Q\  -  *2))  sin  w/2d6,  (B14) 
(Q2  -  Qt)sinhQ2/,e3  +  (*4  -  Q\)s\n*lxch  +  (Q$  -  <?4) sinh Q2/2d3 

=  (jt4  -  Qx)cos7rl2dc. 

The  first  and  third  equations  of  (B14)  yield 

c5  =  d6cotirl2,  (B15) 

and 

c3sinhQ2/]  +  d3sinhQ2/2  =  0. 

The  second  equation  of  (B14)  yields 


,  _  .  ,  _  .  7rd6(l  +  \/3») 

c3cosh<?2/,  -c3coshQ2/2  =  - ■  ■  ■  ■ — . — 

2  Q2  sin  itl2 


Hence, 


__  deir(l  +  v/3i)sinhQ2^2 
2Q2sinW2sinhQ2  ’ 


(B 16) 


dQjr(l  f  V^i)sinhQ2^i 

flg  — - - - * - 

2Q2sinhQ2sin7r/2 


(£17) 


From  the  first  and  second  equations  of  (B13),  we  find  that  Cj  is  the  complex  conjugate  of 
C3  and  that  d\  is  the  complex  conjugate  of  d$. 


APPENDIX  C:  Evaluation  of  inner  products 


We  first  calculate  the  boundary  integrals  that  arise  in  (20)  and  then  the  entries  of 
the  matrix  defined  in  (18).  We  denote  b{  =  (0*,«‘,u >",A*)  and  x°  =  (0 ,u,w,h). 
In  addition  to  the  boundary  terms  arising  in  (20),  we  also  integrate  the  term  arising 
from  p  in  (6,,  L\x®)  by  parts.  This  yields  another  integral  over  the  interface,  which  we 
combine  with  those  from  (20)  into  an  expression  r  tJ  .  The  form  of  V  a  can  be  read  off  from 
the  calculation  of  the  adjoint  in  Appendix  A.  The  terms  remaining  in  (3,,  L\x®)  will  be 
denoted  by  ((6,,  L\X°)).  We  thus  have 

{bitLox})  +  <6,-,L,xJ)  =  (L'0biyx))  +  <(fc,L,xJ))  +  Tij 


where 

_  30  30*  du  dw , 

r‘» =  y,  ®  f al  ■  dT  +  Ft  +  a^> 

IT  ^ 

-©*(-2  Pih~  +  h(Mi  +  —PS))dx 
oz  2 

Here,  the  interval  of  integration  I  extends  over  one  wavelength  in  x,  at  z 
equations  (19)  become: 

<6i,*i>  =  rii  +  «3„L,x°», 

(bi,xl)  =  T|2  +  ({biiLiXt)), 


(Cl) 
/ 1 .  Hence, 


(62, x})  =  1*21  +  ((62,Lix°))  +  (61, x}),  (C2) 

{62^2)  =  ^22  +  {{&2>  ^1*2))  + 

where 

3  2 

Tn  =  2\/2(9P^-fcos  jr/j  +  (Xfj  +  ^-P5)  sin  x/i), 


and 


Tia  =  2>/2(Pircos*7|{(-3  +  9— -)f  -  tf»)  +  (Afri  +  —  P$)sinir/|), 

2  2 

r»i  =2j(f^  +o;(*1  +  jPS)).*i, 

x  2  2 

=  2v/2(^(l-llP)fco6ir/1+2P4(c1(g?-^-)2Q1coeh(?,/I+cs((?l-^-)2(?2coBhQ2/,) 

4  K“  &  « 

.  _._3  7leosinli  ,  .  .  v 

+0  +  **)*;»  -H/mnirfi) 

4  am  *i\ 

+(Mi  +  ^-P5)(cisinhQi/i  +  e3sinh<?2/i  +  (1  +  4  )r~c<Mirli))- 

2  m  OJT 

We  will  see  later  that  r22  is  not  required.  Since  L\x\  vanishes  in  both  fluids,  the  first 
equation  reduces  to  (6i,z{)  =  I'll.  Since  x\  =  -Oi, 

—  —(1  +^)(^i,«i)  +€r„  +  0(e2)- 

We  find  that  (6i,ai)  =  0  so  that 

®ii(£»£)  =  ^11  +  0(c2).  (C3) 

Since  x7  =  -a i  -  o2, 

^ i2(«» =  ~(1  +  d)(6i,o2)  +  e(r i2  -I-  ((hi,l»iz2)))  +  0(f2) 

We  set 

a  =  -1  +  ric1/2  +  r2e  +  0(  e3^2) 

so  that 


^i2(£>^)  —  — +  £(-r2(6i»<*2)  +  Tia  +  <(6i,Lix°»)  +  0(f3^2).  (C4) 


We  have 


(6  i»<*a) 


1  2 n  [x  9Pjr2sina  jtz  x2cos2jtz  .  2  , 

- : - /  — — - —  + - - - +  sin  Kzdz 

sinx/j  a  J0  2(x2  +  a2)  a2 


3y^(l  +  />) 
sin  7rli  ’ 


«*•’*•*&> =  /„'  s'*-*”* 

+  f  sin2  irz(^Pir2?li +  ^PK2*i  +  ^Pn2m -^Pn2(f +  /}))  + SPfhir2 cos2  *zdz 
Jll  2  2  2  2 

0./5  Q  0  |7l 

=  +  -/>irsin2jr/i(f  + 1  —  f  —  0  —  j)). 


Next, 


♦2i(«,o)  =  ~(1  +  6)(&2»ai)  +  c(r  2i  +  (6i,x}))  +  0(e2). 


We  note  that 


(62,«i)  =  (62»^0«2)  =  {io^2’°2)  =  (^l»fl2)- 


Hence, 


^2i(«,o)  =  ~V^rt(^i»®2)  +  €(r  2i  +  (*i,z!»  +  0{e*'2). 


We  find 


^22(^1^)  —  (^2.^2)  ~  (“1  +  v^ri)(^2»«2)  +  O(e) 


=  -(62,01)  -  ^1(62,02)  +  C>(c). 


Collecting  the  0(e)- terms  from  the  equation  det¥,;  =  0,  we  find 


Collecting  the  0(£3/3)-terms,  we  find 


_  fia  +  I^i  +  {(bj, L\x%))  rl  (62,02). 

2  2  <*!,«,)  +  2  1  +  (61,02) 


(C 11) 


We  now  need  to  calculate  {b?,  a^)  given  by 


fl,  r  1  o 


.  2jt  .  [  1  f  .  0‘sinxs  txcosirzu"  ,  . 

\»2»«2/  =  — ; - rt  /  +  /  / 7—5— — 5t  +  - +  tv  sin irz  dz, 

asinir/i  Vo  Jit  (»2  +  «2)  a 


where  63  =  (0*,ti*,tv*,A*)  and  02  is  given  by  (B2).  We  express  the  integrand  in  terms  of 
tv*  by  using  equations  (B4):  t'fi*  =  ^  and  0*  =  3  sin  irz  +  *£  A2  tv*.  Hence, 


2\/2  1  f1  4  p 

(62,02)  =  ^7— |  {^3  +  (y  +J  A2  tv*  sin  xz  +  3tv*  sintrzdz). 


Integration  by  parts  simplifies  the  first  integral  as  follows: 


^  A2tv*sinirzdz  =  sinx/j  [[^“T"]]  +  J  tv*sintrzdz. 


so  that  we  are  essentially  left  with  having  to  evaluate  [[^§-]]  and  /  tZ)*  sin  xzdz  and 
substituting  them  into 


,.  .  2\/2  .  1  4Psin7r/j  rrd3tv 

(621  a2>  =  +  ,fr 


sin  it/ 1  jt2  3tr4 


rrO^tV*  v,  f 

[[-^3-]] +3(1  + P)  J  w’sinnzdz}.  (C12) 


The  former  is  facilitated  by  multiplying  equation  (B5)  by  sin  irz  and  integrating  by  parts. 
This  leads  to 

a»2/ 1  j _ , 

p> 


11  ds3  JJ  4  sin  tr/i 


(C13) 


From  Appendix  B,  we  have 


(1  +  -) 

tv*  =  ci  sinhQis  +  C3sinhQ22  +  cgsinxe - — 


6  ir 


z  cos  irz 


in  fluid  1  and 


t»*  =  disinhQifz  -  1)  +  d3sinhQ2(*  -  1)  +  ^-(1  +  4)U  ~  z)cos?rz 

ok  r 


in  fluid  2,  with  the  coefficients  given  by  (B15)  -  (B17).  After  some  algebra,  we  find 


f1  (1 4-  —  i  1 

/  tD*sin?rzdz  = - — -^(ficotWa  +  — ). 

Jq  12k  2k 


Using  (C12)  -  (C14), 


(62,02)  =  - 


2V2  tP  ,  (1  +  P)2, 


jrsinw/i  Jr  AP 


+  V~  \i~  {hcotxl2+ 
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